Superfluidity and dimerization in a multilayered system of fermionic polar molecules 
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We consider a layered system of fermionic molecules with permanent dipole moments aligned 
perpendicular to the layers by an external field. The dipole interactions between fermions in adjacent 
layers are attractive and induce inter-layer pairing. Due to competition for pairing among adjacent 
layers, the mean-field ground state of the layered system is a dimerized superfluid, with pairing 
only between every-other layer. We construct an effective Ising-XY lattice model that describes 
the interplay between dimerization and superfluid phase fluctuations. In addition to the dimerized 
superfluid ground state, and high-temperature normal state, at intermediate temperature, we find 
an unusual dimerized "pseudogap" state with only short-range phase coherence. We propose light 
scattering experiments to detect dimerization. 
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The long-range and anisotropic nature of dipole-dipole 
interactions offers new opportunities for ultra-cold polar 
molecules, beyond what is possible for cold-atom systems 
with only short-range, isotropic contact interactions []|. 
A variety of exotic many-body states including p^ + ipy 
fermionic superfluids [2|, and nematic non- Fermi liquids 
[3|, arc predicted to occur in cold dipolar systems. Ad- 
ditionally, polar molecules could provide a robust tool- 
box for engineering novel lattice-spin Hamiltonians [J] or 
hybrid devices for quantum information processing [5|. 
Recent progress towards trapping and cooling atoms and 
molecules with permanent electric or magnetic dipole mo- 
ments have opened the door to exploring these exotic 
states of matter experimentally |6| . In order to prevent 
the system from collapsing due to the attractive head-to- 
tail part of the dipolar interaction [7| , it has been pro- 
posed ([8|,[9|) to create stacks of dipolar particles con- 
fined to a set of parallel planes. 

In this Letter, we consider a stack of two-dimensional 
layers of polar fermions whose dipole moments, D, are 
aligned along the stacking direction (z-axis) by an ex- 
ternal field (see Fig. [T]). The dipole-interaction, Vd = 

-pr (1^3^), is purely repulsive between fermions in 

the same layer, and partially attractive (for r < ySz) 
between fermions in different layers. The attractive in- 
terlayer component of the dipole interaction induces BCS 
pairing between layers with adjacent layers competing for 
pairing. We demonstrate that competition between adja- 
cent layers favors dimerization, with pairing only between 
even or odd pairs of layers (Fig. [T|). 

We find three distinct phases: a high temperature dis- 
ordered phase, a fully ordered phase characterized by 
dimerized pairing amplitude and quasi-long range or- 
dered (QLRO) pairing phase in each layer (Fig. 1), and a 
dimerized "pseudogap" phase with only short-range su- 
perfluid correlations. The latter phase is particularly in- 
teresting, since it can only be characterized by a compos- 




FIG. 1. Scliematic representation of competition for pairing 
among adjacent pairs of layers, including depiction of the op- 
tical confinement beam which creates the stack of 2D sheets 
(left), and illustration of one of the two equivalent dimerized 
pairing ground states for a many layered system (right). The 
wavy lines illustrate the proposed light-scattering detection 
scheme discussed below in the text. 



ite four fermion dimerization order parameter. There- 
fore, this phase does not admit a mean field (Hatree- 
Fock) description. This is analogous to spin nematics 
[lO| and charge 4e superconductors 11|, which are both 
phases of strongly interacting fermions that can only be 
characterized by composite order parameters. 

Fermionic Pairing in a Layered System - The action 
for an A^-layer system in terms of fermionic fields ip is 



TV 



^ = EE<k(5r + ek-M)^.,k 



z=l k 

N 



- E E ^Ik'^I'.q-k'^tll'V.Vq-k^^.k (1) 
z,z' — l k,k',q 

where z and z' are (integer) layer labels, ip]^ ^(r) creates a 
fermion with in-plane momentum k and imaginary time 
T in layer z. (The r labels have been suppressed above.) 
Vq'^ is the dipolar interaction between layers z and z', 



Fourier transformed with respect to the in-plane separa- 
tion, for example: 1^ = —D^qe~'^'^. 

By solving the BCS gap equation: A^^k = 

~ X^k' ^-k' (V'2-i-i,-k''/'2,k'), we find that the attrac- 
tive interlayer interactions induce fermionic pairing be- 
tween adjacent layers z and z ± 1 (see supplement for 
details). Interaction between next-nearest layers and be- 
yond is small, and will be neglected throughout most 
of this Letter. To decouple the four-fermion interaction 
term, wc introduce Hubbard- Stratonovich (H-S) fields 
A2(rj_) associated with the pairing order parameters 
(where r^ is the in-plane coordinate), and integrate out 
the fermionic degrees of freedom|12[. Expanding the re- 
sulting fermionic determinant to quartic order (valid in 
the vicinity of the phase transition where |A| is small) 
yields the following Ginzburg-Landau (GL) free energy: 
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J2 [d'r(^K\V^A,\' + r\A,f 
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+u\A,\^ + 2u\A,f\A,,+i\'] (2) 
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where Vl denotes the gradient restricted to the xy-plane. 
The GL coefficients are given by k = glprfft, r — vt, 
and u ~ ^j^ where Q is the Riemann zeta- function, 
t — {T ~ Tc)/T is the reduced temperature, ep is the 
Fermi-energy, and v is the two-dimensional density of 
states (for details we refer the reader to the supplement). 
An important feature of this free energy is that the 
H-S expansion does not generate |9^Ap terms, but only 
terms of the form ^^jAp. The absence of jSzAp terms is 
not an artifact of the H-S expansion; rather, it is guaran- 
teed by particle number conservation for each layer indi- 
vidually. Particle conservation for each layer stems from 
the absence of interlayer tunneling, and formally corre- 
sponds to A^Layers independent U(l) phase rotation sym- 



metries, -02 



./2 



02, of fermion fields 0^ in layer z. 



In contrast to other quasi-two-dimensional systems, such 
as superconducting thin films where the behavior of the 
system tends towards three-dimensional as the film thick- 
ness is increased, two-dimensional Berezinskii-Kosterlitz- 
Thouless (BKT) physics remains important even for a 
large number of layers. 

Mean-Field Ground State - The 2u|Aip|Ai+ip term 
in (2) indicates that adjacent pairs of layers compete with 
each other for pairing. For Mayers > 3, the mean-field 
theory predicts that it is energetically favorable for the 
system to spontaneously dimerize into one of two equiva- 
lent configuration where A vanishes between every-other 
layer: |Aj| = ^ [l ± (-1)-''] Aq (see Fig.l). The situation 
for iViayors = 3 is more subtle, and we defer its discussion. 

Effective Lattice Model for Many-Layer System - The 
above mean-field analysis suggests that the relevant de- 
grees of freedom for a many-layer dipolar system are 
Ising-like dimerization between even or odd layers, and 
two-dimensional XY-like phase fluctuations of the in- 
terlayer pairing order-parameters. In order to describe 
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FIG. 2. Schematic depiction of fully dimerized phase (a), 
in-plane Ising domain- wall (DW) (b), z-axis Ising DW with 
pairing between two adjacent pairs of layers (c), and z-axis 
Ising DW with no-pairing for two adjacent pairs of layers (d). 
Green shading between layers indicates pairing (color online). 

phase transitions in this system, we course-grain the GL 
theory (in-plane) over length-scales below the GL coher- 
ence length ^GL = (K/|r|)^/^, and obtain the following 
effective lattice model [18[ 

^ J2 Ji'^z,^,a,,,) [cos (02,, - 02,j) - 1] } (3) 

of Ising variables ai,z S {il} coupled to XY phase- 
variables 6*2, i ::= ArgA2(ri) G [0, 27r] where z labels phys- 
ical layers, i labels lattice sites in the xy-plane, and 
J{o-z,i,crz,j) = Jo (1 + o-z,i) (1 + (^z.j) /4. 

In the lattice model, (Tz = +1 {az = —1) indicates that 
layers z and z-\-l are paired (un-paired respectively). The 
uniformly dimerized ground state of the multilayer sys- 
tem corresponds to anti-ferromagnetic Ising order along 
the z-axis and ferromagnetic order within the xy-plane. 
Ising domain walls correspond to regions where pairing 
switches between the two equivalent dimerization config- 
urations over a distance of the order of the GL coherence 
length, either along the z-axis or within the xy-plane. 
The coupling constants Ar2 and K± reflect the energy 
cost of deforming the magnitude of the pairing order pa- 
rameter, |A|, to form a domain wall along the z-axis or 
in the xy-plane respectively (see Fig. [2]). 

The coupling J(fTz.i,(T2.j) corresponds to the average 
supcrfluid stiffness p ~ ^jAp in the vicinity of lattice 
site {i, z) and determines the energy cost of twisting the 
phase of the order parameter, 02, i, between sites i and j in 
the same plane. The local stiffness is non-zero wherever 
CT2,i = +1, and zero otherwise J12f. 

The lattice model couplings (Ar2,Ar^, Jq) can be es- 
timated from the GL model. An in-plane dimeriza- 
tion domain wall along the x-direction (Fig. ^}p) cor- 
responds to pairing conflgurations of the form A2(a;) = 



An 



[1 -\- {-iya{x)] where Ag 



|r| 



, and a{x) is a func- 



tion that changes from —1 to -1-1 around x = 0, and tends 



J/T 




FIG. 3. The lattice model phase-diagram, calculated using 
the temperature dependence of the GL coefhcients in (3) and 
plotted in terms of dipole interaction strength D^ /d?' and tem- 
perature T, each measured in units of ef (main figure). Phase 
diagram predicted by the effective lattice model for generic 
model parameters, with Kz/K± = 2. Double line indicates 
first-order transition (inset). 

to a constant away from x = 0. Minimization of the free 
energy with respect to a{x) yields a{x) = tanh (2x/^dw) 

where £dw = a/^y^|13- The corresponding free energy 
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cost per unit length is U^ 

There are two possible z-axis domain wall configura- 
tions, shown in Fig. [2]:,d. To determine their free energy 
cost, we consider a system with periodic boundary condi- 
tions along the z axis, and compare the free energy of the 
ground state to that of the domain wall configurations. 
This yields an energy cost per unit area U)^l^ = ^^ for 
both types of domain walls. Setting the lattice spacing 
equal to ^DW^ the energetics of in-plane and z-axis do- 
main walls are reproduced by Kz = 2K± = -^|p-. In or- 
der to determine the lattice phase stiffness Jq, we equate 
the cost of an infinitesimal phase twist, 9z,j = Oz.i + 59, 
in a fully paired layer {az = 1) to the corresponding cost 
in the GL free energy (Eq. [2|). This gives Jo = ^^. 

The lattice model (Eq. [3]) describes three- 

dimensional Ising spins coupled to many independent 
two-dimensional XY layers. For temperatures near or 
below the Ising transition temperature, the Ising vari- 
ables have large correlation lengths and hence see an av- 
erage over many independent layers of XY spins. With 
this self-averaging property in mind, we decouple the XY 
and Ising variables in a mean-field factorization 

Fa = Kz ^ dziCFz'i - K^ ^ (TziCFzj - h'^CTzi 

{zz')A z,{ij) z,i 

i^XY = - ^ ^ [1 + (-l)Vo]' COs{ez,^ - 9z,,) 



(y) 



where k\''^^ =K^ + ^ (^) and A,B = 



(cos 



,(e/o) 



of the cosine terms in even and odd layers respectively, 
ao = {a)F^, and h={^)^. 

The decoupled Ising model and XY models can then 
be analyzed separately but self-consistently. A mean- 
field analysis is adequate for 3D Ising model. The phase 
action is treated by a variational self-consistent harmonic 
approximation (SCHA) [13|. While the SCHA provides a 
reasonable estimate of the location of the 2D BKT transi- 
tion, it spuriously predicts a strong first-order transition 
in which (Fxy)scha drops abruptly to zero at the XY 
transition temperature, Txy- At higher temperatures, 
the SCHA dramatically underestimates the contribution 
to energy density from phase fluctuations. In order to 
avoid this undesirable feature, we supplement the SCHA 
value for (cos Ay6')scHA with a high-temperature expan- 
sion for T > Txy- 



{cosAijO'-^^) 



(cosAy6i(^ 



/SCHA 



J(ao,ao)/2T 



;T<TxY 
■,T>TxY 



(4) 



Of ))fxy ^ 1 ^I'C the averages (with respect to Fxy) 



Fig. [3] shows the phase diagram predicted by the ef- 
fective lattice model. The main figure displays the phase 
diagram where the model parameters are taken from the 
GL coefficients in (3). The BCS transition temperature, 
T^^^, is obtained by solving numerically the BCS gap 
equation for the dipole potential. Whereas the dimeriza- 
tion transition occurs close to the mean-field BCS tran- 
sition temperature, T^'~^^, the BKT transition to phase 
QLRO occurs at a lower temperature, leaving an interme- 
diate region with full dimerization but only short range 
superfluid correlations. 

Recent experiments on 3D clouds of ultra-cold 
^°K® Kb molecules have achieved densities on the order 
of n^d = lO^^cm"'^ and permanent electrical dipole mo- 
ments of up to 0.566 Dcbye [6i]. If similar densities were 
achieved in a layered system with layer spacing on the 
order of 400nm, the ratio of typical dipole interactions to 
Fermi-energy would be D'^ /{Aireod^ep) ~ 3. 

While the GL parameters in Eq. [2] provide an ini- 
tial estimate of the lattice-model coupling constants, in 
principle, the model coefficients can be renormalized by 
higher order terms in the GL expansion. The inset shows 
the phase diagram for generic values of the model param- 
eters K± and Jq with Kz/K± = 2 (the qualitative fea- 
tures do not depend sensitively this ratio). An additional 
feature emerges for generic coefficients: for J sufficiently 
bigger than K, there is a tri-critical point where the BKT 
and Ising transitions fuse into a weakly first-order phase 
transition. 

Order parameter and detection - The dimerized phase 
breaks translational symmetry in the z direction. It 
can be characterized by the following four fermion or- 
der parameter: V = {nz-i^rnz,r — '^2,r"-z-i-i,r), where 
nz,r = "01 r'^z,r IS the local fermion density. For finite 
transverse confinement, in the dimerized phase, every 
two paired layers shift slightly towards each other. The 



displacement scales as bz ex Oj^, where fi^ is the layer- 
confinement frequency in the z-direction. The dimerized 
phase can be detected by the appearance of new Bragg 
peaks in elastic light scattering (see Fig. [1]) with wavevec- 
tor Q = nTTz/d, n = 1,3,..., with intensity ~ (5z^. 

In the strong-confinement limit, £7^ — ^ c», the particle 
density does not show any sign of dimerization. How- 
ever in this regime, the dimerized phase could still be de- 
tected by measuring correlations between the amplitudes 
of light scattered at different wavevectors: (nqTiq/) ex 



"O^q-Hq' 



5Q_q_q'I?, where q and q' are two scattering 



wavevectors and no is a constant. 

Three Layer Case - The three layer system is a spe- 
cial case that requires more careful analysis. If one pro- 
ceeds as above and includes interactions only between 
neighboring layers, the system possesses an extra SU{2) 

symmetry generated by: I^ = J (Pr (ipJ^tp^ — iplipi 
and /^ = J d'^r (iplipi^iTplTpsj. The [/(I) generator 

N2 = J (frij;l%jj2 completes the SU{2) symmetry to U{2). 
These generators commute with H = ?^kin + V12 -f V23, 
where Vij is the interaction between layers i and j. 
This U{2) symmetry dictates that, to all orders in the 
GL expansion, the free energy should be a function of 
(IA1P + IA2P) only, which does not energetically distin- 
guish dimerization from uniform pairing. 

However intralayer and next-nearest neighbor inter- 
actions V = V13 -I- X]7=i ^jj break the SU(2) symme- 
try of the three layer system, and generate a quartic 
term of the form — |w|]A_ip|A2p in the GL free energy. 



This term is relevant [14| (in the rcnormalization group 
sense) and hence, we expect the trilayer system to ex- 
hibit uniform pairing with |Ai| = IA2I. In contrast, for 
-^Layers > 3, already the dominant nearest neighbor inter- 
actions strongly favor dimerization and V only produce 
small sublcading corrections. 

Discussion - Our analysis of the layered dipolar Fermi 
system predicts a sequence of two phase-transitions: an 
Ising-like dimerization transition followed by a BKT 
transition to phase QLRO. One can generalize to one- 
dimension, and consider a stack of one-dimensional tubes 
of dipolar fermions [15| . In this case, no phase ordering 
can occur at any finite temperature, since the phase dy- 
namics are strictly one-dimensional. However, a dimer- 
ization transition is still possible, leading to wider range 
of dimerized, non-supcrfiuid phase[16|. 

We expect that the Ising-XY model description of the 
layered dipolar fermions will be insufficient deep in the 
BEG regime where interaction energies are dominant 
compared to the Fermi-energy. For sufficiently strong 
interactions or sufficiently dense systems, the system will 
form a Wigner crystal [l7[. Another possibility is that 
the formation of longer chains of three or more dipoles 
may become important [9|. In a regime where chains 
of n dipoles are favorable, a many-layered system would 



undergo n-merization rather than dimerization. Gorre- 
spondingly, an n-merized phase may undergo an n-state 
clock-model-type phase transition which generalizes the 
Ising-type dimerization transition considered above. Fur- 
thermore, for even n, bosonic chains could condense into 
an exotic superfluid of dipolar chains. Such states offer 
an intriguing chance to examine the relatively unexplored 
boundary between few-body interactions and many-body 
phase transitions, and deserve further study. 
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Appendix A. BCS Gap Equation for Bilayer 

In this section, we show that the attractive component 
of intralayer dipolar interaction induces BCS pairing in 
a bilayer system. We ignore repulsive intralayer interac- 
tions, as these serve only to renormalize the Fermi liquid 
parameters of the dipolar system. The interaction be- 
tween dipoles in adjacent layers separated by interlayer 
spacing d along the z-axis, and by distance r in the xy- 
plane, and with dipole moments D polarized along the 
z-axis by an external field is: 



Vdip{r) = 



D' 



(^2+^2)3/2 \^ ^2 + 



1- 



3^2 



rf2 



(5) 



where we work in units with 47r£o = 1- Fourier trans- 
forming ([5|) with respect to in-plane coordinate r, one 
finds: 



"^9 



-D^qe 



-qd 



(6) 



We assume that pairing in the s-wave channel domi- 
nates, and that the transition temperature is set by con- 
densation of Cooper pairs with zero center of mass mo- 
mentum. The self-consistency equation for the BCS or- 
der parameter A at temperature T reads: 



A, 



X] ^k-k^l^^ (V'z+l,-k''02,k' 




FIG. 4. (color online) Numerical solutions of the BCS gap 
equation (O for the nearest-neighboring layer dipole poten- 
tial ((6]), at fixed particle density, and with d/Xp = 1, where 
\f ~ 211 /kp is the Fermi wavelength, (a) T = BCS pairing 
as a function of interaction strength D^ jS (where energies 
are measured with respect to the Fermi energy ef of the un- 
paired system). Solid line shows solution for A, and dotted 
line shows chemical potential [i,. The BCS regime (/x > 0) and 
BEC regime (/x < 0) are highlighted in red and blue respec- 
tively, (b) Temperature profile of the BCS gap for interaction 
strengths: D'^ /d'^ep = 2, 1.75, 1.5, 1.25, and 1 (highest to 
lowest respectively). 



In addition, one of us (D.W. Wang) has performed 
detailed studies of BCS pairing in A/" = 2, 3, and 4 layer 
systems. These studies confirm that interlayer pairing 
occurs in these few-layered structures, and demonstrate 
explicitly that dimerization is favored for A/" = 4, but 
not for M = 3, in agreement with our analytic results 
from the Ginzburg-Landau free energy. This work will 
be published elsewhere. 



iV^i^izlll_^tanhr^ 
2^^ Ek, \ 2T 



(7) 



In (HI), £;fc = ^Cfc + A^k where ^k = ^ - ^^, m is the 
effective mass, and fi is the chemical potential. Since 
the number of particles on each layer is fixed, one must 
simultaneously solve for A using ([7|) and for fi by fixing 
the particle density n: 



E 



l_i^tanhf:^ 



(8) 



We solve (|7l8p numerically, using the expression ^ for 
the dipole potential, and find a non-zero solution for any 
dipole moment D. Fig. (|4])a. shows the T = solution 
for A as a function of interaction strength D^/d^ep, and 
Fig. (|4])b. shows the temperature profile of the BCS gap 
for various interactions strengths. These results demon- 
strate that intralayer interactions induce pairing for any 
value of dipole interaction strength. 



Appendix B. Derivation of Ginzburg-Landau 

Free-Energy 



In this section, we include for completeness, details of 
the derivation of the Ginzburg-Landau (GL) free-energy 
^ from the microscopic fermionic action ([T]) for an Af 
layer system. This derivation, based on the Hubbard- 
Stratanovich (H-S) transformation, follows a standard 
route to deriving the GL free-energy in BCS theory (see 
for example A. Altland and B. Simons, Condensed Mat- 
ter Field Theory (Cambridge University Press, Cam- 
bridge, 2006) p. 271-309). 

Starting with ^ we introduce H-S fields 

A^(ki,Q_L) = -Ek'^k-ll^|^^(^'^+i,-k'+q/2'/'2,k'+Q/2)- 
The subscript z labels layer, and the two momentum 
labels (kx,Q±) refer to the relative displacement and 
center of mass motion of the Cooper pairs respectively. 
The inclusion of Q^ 7^ allows for spatial varations of 
the order parameters A^. 

With the introduction of the H-S fields A,, the 



fermionic action now reads: 



z,k,k' ,Q 



Ha 



A.,fc,g {V~^)^^^, A.,.',Q + Y. *t^^HA*fc,Q 



Al,A:,Q 





Ai,fe,Q _ 
-A2,fc,Q iw + C 



(9) 



(10) 



7 



where \l/fc,Q = (V'i,-fe+Q ^"2,*: ■03,-fe+Q V'4,/c • ■ •) , and 
Lo is an (imaginary) Matsubara frequency label. Here, 
{V~^)k.k' denotes (fc, k') component of the inverse (in the 
operator sense) of the dipole potential ^ . As will be ex- 
plained below, it is not necessary to explicitly compute 
this inverse in order to develop the GL theory. Since we 

I 



are concerned with finite temperature phase transitions, 
we neglect quantum fluctuations by treating A as inde- 
pendent of uj. This corresponds to developing the GL 
free-energy for the w = component of A^. We also as- 
sume that the relative displacement momentum k profile 
of Afe_Q does not fluctuate substantially from the mean- 
field form. This corresponds to fixing the form of the 
Cooper pair wave-function to the one most energetically 
favored at the mean-field level, and is justified by the fact 
that the dominant instability towards pairing will occur 
with this pairing profile. With this assumption, the k 
labels on A are non-dynamical and will be dropped from 
subsequent expressions. We have validated this approach 
by checking that the resulting GL free-energy derived in 
this way can accurately reproduce the results of the nu- 
merical solution to the BCS gap-equation. 

Integrating out the Fermions gives the following effec- 
tive action for the Hubbard-Stratonovich fields A,: 



5[A] = Y^ 5(2n) j^j ^ rp^lj^ [^ ^ Go ^Ha] 
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r 



where we have dropped all irrelevant constant terms that 
do not depend on A, and Gg^ = ^^, G^^ = j^ 
are the particle and hole Green functions respectively 
(note the precise form of the bottom right entry of the 
second line of (|lip depends on whether M is odd or even). 
Also, in pT|) momentum labels on A and Go have been 
suppressed where possible in order to conserve space. 

We truncate the above series at quartic order, which 
is formally justified near the pairing transition tem- 



perature T^^ where A is small. The GL coefficients 
{k,,t,u\ can then be computed by explicitly evaluat- 
ing the trace over the products of fermion Green func- 
tions appearing in (|lip . The terms quadratic in A 
are: 5^2) = Y.q^Q'\^{Q)? where T{Q,T) = V'^ - 
^^ Go^pGq„^q, and where Vt is the system volume. 
The mass term, with coefficient r, determines the en- 
ergy of having a uniform pairing amplitude |A|. To com- 
pute r, we note that at T^'^ we have: r(0,T^'^'') = = 



V-^- ^ Ep Gg pG[j__p, indicating that for T near T^^ 
d 



n 
one can expand 



(rp _ rpMF\ ]^ \ ^ ^p i^h 
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where np{e,T) is the Fcrmi-distribution at energy e — fi 
and temperature T. The coefficient k, of the gradi- 
ent term KJVxAp, is obtained by expanding r{Q,T) to 
quadratic order in Q: 
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^ -dlnpiip^T) _7C(3)£F 
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Note that obtaining k and r by expanding V{Q,T) to 
leading order in Q and T respectively has allowed us to 
neatly sidestep the explicit computation of V~^ . Anoth- 
ering interesting feature is that the details of the dipolc 
potential are fully contained in the single parameter, 
T^^, which we obtain by numerically solving ([7]). 

Finally, we turn to the evaluation of the quartic term 
coefficient u, which comes from terms of the form 5^'*^ ~ 
|A|^ ^(G^G'ip) . The Matsubara frequency summa- 
tion in ^ {G^fi'!_p) can be done explicitly: 



E. 



1 



1 



(ia;„ - ^p)2 {iojn + £,pY 
-dz P 1 



2^1 eP^ + 1 (z - Cp)^ {z + ^pY 
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-^ tanh(/3^p/2) (14) 



8C2cosh^(/3Cp/2) 4^3 

where /? = \/T is the inverse temperature. The remain- 
ing integral over p can be rendered dimensionless and 
computed numerically yielding: 



K = 



32 



dx 



1 tanhxX 1.7 i^ 

r2 cosh' a; ^ a;3 j " 32 r3 

(15) 
Using the above computations of {K,r, u} it is now a 
straightforward matter to evaluate (jlip , which for infinite 
number of layers gives: 
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The relative factor of 2 between the \^z\^ and 
jA^pjAz+ip terms can be obtained by a careful ac- 
counting of combinatorics. Alternatively, this factor may 
be checked by considering a trilayer system whose mi- 
croscopic SU(2) symmetry (described in the main text) 
guarantees that the quartic terms in the GL free energy 
be of the form: S^^Ls = "(lAip + A2|2)2. 



Appendix C. Coupling Between Domain Walls and 
Phase Fluctuations 



In this section, we discuss some subtleties involved in 
coarse graining the GL free-energy ([2]) to arrive at the 
effective lattice model ([3]). In deriving the value of the 
XY stiffness J {at, aj), which couples the dimerization or- 
der parameter to phase fluctuations, we have made the 
simplifying approximation that the the superfiuid stifness 
drops abruptly to zero at the location of an in-plane Ising 
domain wall. 

In reality, there is some residual superfluid stifi^ness 
that varies continuously across the length of the do- 
main wall. Keeping track of this residual stiffness cor- 
responds to including higher order gradient terms such 
as K^^^V ±Az{r)\'^, which includes terms of the form 
|Ap|Vi.A(r)p|V6'(r)p that couple domain walls and 
phase fluctuations. However, these terms are subleading 
and have a negligible effect on the systems phase dia- 
grams. Formally, this is because such higher order gra- 
dient terms are highly irrelevant in the renormalization 
group sense. 

We have also conducted separate simulations that in- 
clude residual phase stiffness at the dimerization domain 
boundaries, and have explicitly conflrmed that the re- 
sulting phase diagram is highly insensitive to the inclu- 
sion of such terms. Specifically, in addition to the usual 
Ji<7i = +l,o'j = +1) = •^o and J(— 1,— 1) = 0, we al- 
lowed for J(+l,— 1) = J(— 1,+1) = Jros, and verified 
that the choice of value for Jics had little discernible ef- 
fect on either the dimerization or BKT phase transitions. 

A second simplifying assumption is used in parameter- 
izing the in-plane domain wall as A(a;) = ■'-2^(1 + a{x)). 
This parameterization implicitly assumes constant phase 
over the range of the domain wall, which is justified be- 
cause ^DW ^ ?GL. Since ^gl is the shortest scale over 
which the phase is well-defined, it is essentially constant 
on lengthscales < ^gl- 



